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ABSTRACT 

This paper uses the Schwarzschild metric to derive an effective refractive index and 
acceleration vector that account for relativistic deflection of light rays, in an otherwise 
classical kinematic framework. The new refractive index and the known path equation 
are integrated to give accurate results for travel time and deflection angle, respec- 
tively. A new formula for coordinate acceleration is derived which describes the path 
of a massless test particle in the vicinity of a spherically symmetric mass density dis- 
tribution. A standard ray-shooting technique is used to compare the deflection angle 
and time delay predicted by this new formula with the previously calculated values, 
and with standard first order approximations. Finally, the ray shooting method is used 
in theoretical examples of strong and weak lensing, reproducing known observer-plane 
caustic patterns for multiple masses. 

Key words: gravitation - gravitational lensing: micro - methods: numerical - accel- 
eration of particles - planets and satellites: detection - gravitational lensing: strong. 



1 INTRODUCTION 

The deflection of light by massive bodies has been of interest to mathematicians and physicists from time to time since Newton 
suggested the possibility IjNewtonl ([1704])). This deflection was calculated in the late eighteenth and early nineteenth centuries, 
treating light as a classical particle. The deflection was again calculated by Einstein in the early twentieth century, using his 
new general theory of relativity, to be twice the previous classical result. The measurement of the deflection of light passing 
close by the sun was widely publicized as a dramatic confirmation of general relativity, in the now famous 1919 expeditions. 

In the last three decades, gra vitational l ensing has become an important tool for astrophysicists, especially in searching 
for dark matter (first suggested in lPaczvnskil 1 19861 )^ and exoplanets. By 1991, astrophysicists were sug gesting that exopl anets 
could be found using microlensing, and since 2004 at least ten planets have been found in this way (jSumi et 3 (|2010h . In 
microlensing, light from a background star passes close to the lensing system, and is deflected around the lens. Because of 
this, more light rays reach the observer, producing magnification of the background source. This magnification changes over 
time, as the source, lens and observer move into and out of alignment. The details of the magnification over time are plotted 
in a 'light curve', which is simply intensity versus time. A planet orbiting a lensing star can make changes in the magnification 
pattern at the observer's plane. Such changes show up as variations to the shape of the simple light curve. 

Interpretation of these light curves is difficult, as this is an inverse problem, in which observers seek to determine the 
details of the lensing system which gave rise to the observed data. In particular, researchers are often trying to find planets in 
the lensing system, and to determine characteristics of such planets, primarily orbital radius a nd ma ss. Details on reproducing 
a model of such a m ulti-body le ns were outlined at least as early as 1996 ( Wamb sganssl 1)19971) ). For a brief history see 
Schneider et all (119921 ) . pages 1-2. iMao et all (|2007l ). or the recent review by lEllisl (120100 . In this paper, microlensing will be 
considered, although the formulae derived here are just as valid for macrolensing situations. 



1.1 Lensing Models 

It is customary ( Wambsganssl ( 1997 )) to use a 'thin lens' model, in which the effect of the lensing system is confined to 
the plane containing the lensing objects (the 'lens plane'), this plane being normal to the line from observer to lens (or 
alternatively, source to lens). A deflected light ray thus consists approximately of two straight lines, with an abrupt angle 
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Figure 1. Thin Lens Approach and Typical Light Curve 




change in the lens plane. The magnitude of this change is given by a simple formula: AO = 2r s /ro where ro is the point of 
closest approach of the light ray to the star or planet and r s is the Schwarzschild radius, r s — 2MG/c? . Here, M is the mass 
of the star, G is Newto n's gravitational constan t and c is the speed of ligh t in a vacuum. 

Current methods ( Zabel fc Peterson ( 20031 ) and Wambsganssl (1997)) shoot rays from the observer to the lens plane, 
deflect by the angle as calculated above, and then draw the ray from there to the source plane. Equivalently, light rays may 
be followed from the source to the observer, mimicking more closely the actual physics. The density of rays at the source 
plane (alternatively, the observer's plane) is thus mapped. By tracing various linear paths across this 'magnification map', 
to simulate the relative movement of the source star, light curves are generated by plotting density of rays as a function of 
time. Fig. (1) shows schematically the method for mapping the amplification for an observer travelling at constant speed in 
the observer's plane. Simple light curves such as the one in the right of Fig. (0 provide information about the mass of the 
lensing object, provided that the distance to lens and distance to source can be estimated. If the lens is a binary system, there 
are deviations from this simple light curve, which may yield information about the orbital distance and relative masses in the 
system. 

Thus the purpose of a microlensing model is to produce a magnification map due to light being deflected by the lensing 
system. Note that every such magnification map allows for arbitrarily many light curves, depending on the location of the 
observer's (or source's) path across the plane. Because of this, it is, in general, a difficult matter to find a model to fit an 
observed light curve for a lens involving more than one mass (that is, a bina r y or p lanetary system). The problem becomes 
much more difficult when the lens involves more than two objects ( Mao et all 1 2007D ). 

The approach presented in this paper is to generate a magnification map by shooting rays through a lens of smoothly 
varying refractive index. The path will be curved rather than two straight line segments. Even so, the majority of the deflection 
occurs very close to the lens plane. 



1.2 Kinematical Approach 

It is the purpose of this paper to consider gravitational lensing using an alternative approach. This approach can be described 
as an equivalent classical kinematical formulation, which nevertheless replaces Newton's formula for gravitational acceleration 
with a formula derived from general relativity. Using this kinematical approach, each curved light ray is traced from the source 
to the plane of the observer, rather than using the simple angle change formula commonly used in lensing models. 

We will start by considering the path predicted by general relativity for a massless test particle (a photon) near a massive 
body. The changes in angle along this path will be integrated to give the deflection of the light ray predicted by general 
relativity. As a check, we will compare this result with known formula for approximating the deflection to first order in the 
relevant small constant. Next we will derive an 'effective refractive index' due to the massive body. Integration will be used 
to predict the Shapiro delay, which will also be checked against the known first order prediction of general relativity. A new 
acceleration formula will then be derived and tested against these known values for the deflection and delay using a forward 
integration method. Finally, we will describe and exemplify a procedure for using the new acceleration formula to produce 
magnification maps for multi-body lenses. 
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2 LIGHT PATHS IN A SCHWARZSCHILD SYSTEM 

The Schwarzschild metric describes spacetime outside a non-rotating, electrically neutral, spherically symmetric mass density 
distribution M and is taken as a valid approximation for local spacetime structure in the vicinity of any massive body (stellar 
systems, including black holes, but also planets) having negligible charge and angular mome ntum. In spherical coordinates 
the metric components are given in terms of the invariant interval dl as l|Misner et al.l (jl973h ) 

,,,2 _r~ Ts_ c 2 dt 2 r _ dr ? _ r 2( sin 2 8d(j) 2 + M 2^ ^ 



r r — r s 

where r s = 2MG/c 2 is the Schwarzschild radius, G is Newton's constant, and c is the constant speed of light in vacuum. Here 
r , 8, <f> are coordinates which correspond to standard spherical coordinates in the reference frame of an observer at rest far 
from the system. 

Light travels on null geodesies with dl = 0, so equation |T]) becomes 

cdt = — J dr 2 + r(r- r s ) (sin 2 9dtj> 2 + d9 2 ). (2) 
r — r a V 

Applying Fermat's principle, that light follows a path that extremizes travel time T, we can consider the functional: 



dt 



- f —^—Jdr 2 +r(r- r a )(sm 2 6d<t> 2 +d8 2 ). (3) 
c J r — r s V 

To illustrate the use of equation © we firstly derive a path equation which is equivalent to the usual trajectory equations 
for null geodesies, and numerically integrate it to give the deflection of a light ray according to general relativity. 

Without loss of generality, for a single light ray, the coordinates can be oriented so that the ray is in the plane 8 = n/2. 
Then dd — in equation ((3), which consequently may be re-arranged to give 



c J r — r 

adopting the notation <j>' = d<j>/dr . Let F be the integrand in equation that is 



sjl + r{r - r a )<j) l2 dr (4) 



F(r, <f>, </>') = Vl + r(r-r s )^' 2 . (5) 

Then the Euler-Lagrange equation is 

dF d OF 

dj-d-r { W ) = °- (6) 

It can be seen in equation ©, that there is no explicit dependence of F upon <f) so clearly dF/d(f> = 0. Therefore this term 
vanishes in equation (0 which then has the immediate first integral 

dF _ r 2 4>' _ 

d(j}' ~ y/T+Hr -r s )(f)' 2 ~ ' 
in which K is & constant. Now l/<f>' — dr/d(j>, so rearrangement yields the first order separable ODE 

dr 



dep -V KS-< r - r ->- 

This constant K can be determined. At the point of closest approach to the mass (call this point r = ro), the radius is at a 
minimum, that is, dr/d<f> = 0, so it follows that: 

K 2 = ^. 

ro - r s 

Thus the path is defined by: 



dr r 4 (r - r s ) 

dt = v — ^ — 

It can be easily seen th at substituting u = 1/r followed by differentiation gives the well known second order equation as 
1 Capozziello et all jl99it l'l 

d 2 u 3r s 2 , Q s 

From equation (O, an integral can easily be written to evaluate the total change in of a light ray passing a massive 
object, which will be twice the change from the perihelion out to infinity: 



A<f> = 2 



dr 



\/r 4 (ro - r a )/rl - r(r - r a ) ' 
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This elliptic integral cannot be evaluated with a finite number of simple algebraic terms. Numerical integration could be used 
to evaluate the deflection. However, the integrand is infinite at r — ro, making the accuracy of any numerical evaluation 
questionable. By using a substitution r = 1/ cos(2?/;) , the singularity can be removed. After simplification we obtain: 

ir/4 

' /2ro f # 




\J 2ro /r s +2 — 4 cos 2 tp — sec 2 ijj 

This integrand is well behaved over the interval, so it can be numerically integrated to any desired accuracy, for example, by 
gaussian quadrature. An undeflected ray will have A<f> = ir, so the deflection for a ray will be S = A<f> — n. The deflection is 
usually very small, so in order to avoid a 'loss of significance error', this subtraction should be performed in the integrand. 
For readability, let a = 2ro/r a , and let /? = 2 — 4 cos 2 ip — sec 2 ip. In this notation, the deflection 8 can therefore be expressed 
as 

7T/4 

8 = 4^ (^J= --!=)#• 

J sja + p y/OL 



The two terms in the integrand are combined, so that the subtraction can be performed explicitly. This yields: 

vr/4 

5 = 4 / ( ~P (9) 

o 

For the purposes of this paper, we take the solar radius to be 696000 kilometres, and the solar Schwarzschild radius to be 2.95 
kilometres. Numerical integration of equation ([9]) for a ray passing near the surface of the sun (r s /ro — 2.95/696000))gives a 
deflection of 1.74851634161261 arcseconds. As the path equation contains all of the information about the general relativistic 
path of the photon, this is the deflection predicted by general relativity to the level of precision shown. This deflection angle 
will be use d later to confirm the accuracy of the kinematic approach. As a check, we can consider Einstein's estimate for the 
deflection ( Schneider et all 1 19921 ). page 3) 



2r 

—L = 1.74850913341648 arcseconds. (10) 
ro 



This estimate is found to correspond to our calculated value to first order in r 3 /ro, as expected. 



3 NEW REFRACTIVE INDEX AND TRAVEL TIME DELAY 

The approach presented here makes use of an expression for the refractive index, n. The Schwarzschild metric from section 
2 will be used to derive an 'effective' refractive index. As suggested above, the functional $Jfy can be arranged in the form 
T — J dt = J(dt/ds)ds, where s is an arbitrary parametrization of the ray path. By choosing ds to be the element of arc- 
length along the path, the speed is then v — ds/dt, and thus the refractive index is n = c/v — cdt/ds, and we obtain finally 
T = 1/c J nds, with the refractive index then having the form 



■\Jr' 2 + r(r - r s )(sin 2 6(j>' 2 + 6' 2 ), 



where r' = dr/ds, 9' — d9/ds and 4>' = d<f>/ds. (Note that there is no suggestion here of any physical effect on the speed of 
light - indeed, any local measurement of coordinate velocity is guaranteed to result in the usual speed c). 
Taking again the two dimensional case, with 9 = tt/2, the refractive index has the form 

n=— — Jr 12 + r(r-r s )(j>' 2 . (11) 
r — r a 

As an example, the delay can be calculated for light to travel from an object at earth radius, sk im past the sun , and back 



out to earth radius. This is the calculation needed for the radar echo delay test of general relativity ( Shapiro! (1964)). (In fact, 
for that particular test, it would be necessary to consider another planet such as Venus or Mars, and the calculation would 
need to be performed for each leg of the journey. For the purposes of this paper, it is simpler to imagine a reflecting satellite 
at the same orbit as the earth). The problem can be pictured as in Fig. [2] not to scale. The deflection is small, so that the 
path appears as a straight line. In fact the curved path that the photon takes is derived above in equation (0. The path is 
symmetric about the point of closest approach (ro), so that exactly half the delay can be obtained by integrating from ro to 
r e , the radius of the earth's orbit. If the sun had no effect on the light path, the distance from perihelion to earth would be 
the straight line distance \Jr 2 — r\. 

The delay can be calculated using the new refractive index, equation (|11[) , and multiplying by ds so that 



nds = — \J dr 2 + r (r — r s )d4> 2 = — \j 1 + r(r — r s )(^-) 2 dr. 
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Again, (d(f>/dr) 2 is given by the path equation (J7J. The time taken for the trip will be twice the time to go from the point of 
closest approach to the sun, ro to the earth's orbit, r e , that is 



c 



r Jr 



t = i r n d s = - r ^-Ji + ^ r{ :~ rs) , ^ 



cj r~r a V r 4 (r ~ r s )/rg - r(r - r s ) 



The integrand is infinite at r — ro, so before integrating, we make the substitution cos(V>) = ro/r. After rearrangement we 
get: 



T = 2r§ Co/^) # / (r -r.)(l + co8V>) 

c 7 cos 2 i/>( r o — r s cos VO y (^o — r s )(l + cos?/)) — r s cos 2 tp ' 

This integrand is perfectly well behaved over the interval, so it can be integrated to arbitrary precision by using, for example, 
Gaussian Quadrature. As for the deflection, the delay is very small, so it is important to subtract the straight-line time before 
integrating. The straight-line time can be written as: 

arccos(rQ /r e ) 



2r f^cos(r 0/ r s) ^ 



t °=cV^=t7 cw ^ 



The delay is obtained by subtracting equation (fT3l) from equation (fl2ll to give: 



At 



2ro 
c 



/" 

Jo 



s(ro/re) (1 - y)d^ 
7 cos 2 tp 



(14) 



where 



7= (l- (r s /r ) cos i>)\ 1 



r s cos 2 tp 



(ro — r s )(l + cost/;) 

This expression (|14|) is now in a form that minimizes errors caused by subtraction of large terms. Performing this calculation 
in Matlab gives a delay of 129.0896086 fis. As for the deflection, this calculation can be performed to arbitrary precision, and 
accurately describes the delay predict ed by general relativity. 

For comparison, Weinberg 1 1972h . pages 201 - 203, gives the following formula for the delay, to first order in r s /ro: 



2 In 



r e + \Jrl+ r% 



ro 



ro 



r e + r 



(15) 



This gives a delay of 129.0894053 /is. The fractional difference between this and the delay calculated above using equation 
(fT4l> is: 



129.0896086 - 129.0894053 



1.57x10" 



129.0896086 

which is of the same order a s r s /rp = 2.95/6 96000 ~ 4 .24xl0~ 6 , as e xpected. Note that equation (|15p is usually further 
approximated and simplified ( Weinberg! ( 19721 ). page 203, Kenvon (1990), pages 95-96) by saying that ro << r e . This gives: 



At 



\ r 

which is significantly less accurate, giving a delay of 129.1350325 [is. 
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We now turn to a kinematic approach, and consider differential equations that relate position, velocity and acceleration. 
In a Newtonian system, the acceleration would be given by Newton's law of gravitational attraction, g = —c 2 r s /2r 2 e r . This 
simple formula is not appropriate for the present application, and a new form will now be derived from the present relativistic 
approach. 



4 A NEW ACCELERATION FORMULA 

By combining the metric equation |T]) with the path equation {7} for a photon in a Schwarzschild orbit, the velocity and 
acceleration of the photon due to the nearby mass are now derived. Here, the meaning is that of a coordinate acceleration. 
As a freely falling particle, the photon does not experience any locally measurable force. Beginning with the path equation 
([7]), and making the substitutions 



A 



ro 



and fj,(r) 



(16) 



the path equation becomes 



dr 



Ar z 



(17) 



Next, considering equation |T}, setting dl 2 = 0, 6 = ty/2 and dividing both sides by dt 2 , we have 

cV = r 2 + mW>) 2 (18) 

where r = v r is the radial velocity component and r(j> = v$ is the tangential velocity component. Using the path equation 
(I17[) . we can solve for v r and v$ in turn, to get: 



V4, = r(j> — ±- 



CfJ, 

I Ar 2 



(19) 
(20) 



Thus the velocity vector of the photon along its path is 

v — v r e r + v<f,e<p — ±cfi 



l--/^e r ± 
Ar 2 



1 



'Ar 2 



To determine the acceleration vector, take the derivative with respect to time: 

a — v r e T + v r e r + v'^e^, + v^e^ 

In polar coordinates, the derivatives of the unit vectors are eV = 4> e -<l> an d e </> — ^4' e r, and so the acceleration components 
a r and are 

a r = v' r — (pv,),, 
= v<f, + 4>v r . 

Differentiating v T in equation (|19[1 yields: 



1 ^_ , M(2r - 3r s ) 

Ar 2 2Ar 2 y / l- fi/Ar 2 



Equations (|19|l and ([20} are now used to eliminate the square root terms, so that v r simplifies to 



r s v 2 , 3 2 



The radial acceleration component is therefore 

Or = Vr — (ftVtj, — v r 



V l-+L 

r 2 



2 o 

r s v r 3. 2 
— + (r--r s K 



The acceleration must be related directly to the Schwarzschild radius r s of the mass. There are two terms in equation (|21[) 
that do not have an r s coefficient. These two terms cancel if the positive sign is chosen. Thus, the correct form for the radial 
acceleration is: 



3v d 



/' 

A similar treatment for tangential acceleration component yields 



£10 = Vcf, + <fw r 
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Approach 


Deflection (arcscc) 


Delay (/is) 


Gen. Rel. prediction by integration along path (cans. [51 and 1 1411 


1.74851634161261 


129.0896086 


Usual approximations to first order in r s /ro (eqns . flOl and ITSt 


1.74850913341648 


129.0894053 


Forward integration using new acceleration formula 


1.74851634 ± 10" 8 


129.089609 ± 10" 6 



Table 1. A Comparison of the new method with accurate predictions from general relativity, and the common first order approximations. 



Thus the acceleration vector for a photon near a Schwarzschild mass is: 



-11 



3». 



4, 1 , V r Vs 

_*] Cr + __e* 



(22) 



It is interesting to note that the radius for light to remain in a circular orbit about the mass can immediately be derived 
from this acceleration vector. In a circular orbit, there is no radial velocity, and so 

Sv 
~~2 



(23) 



In addition, an object moving in a circular orbit in a classical kinematical framework has centripetal acceleration vector 



When equations (|23|l and (|24[1 are equated, we obtain 



3r s 



(24) 



(25) 



Thus, a photon at the '3/2' radius given in equ ation (|25l) is tr apped in a circular orbit about the mass, in accordance with 
the known result predicted by general relativity (|Carrolll (|2004r ). p 212). 



4.1 Kinematic Ray Shooting 

With known acceleration components, it is now possible to set up a standard system of differential equations for ray tracing 
in polar coordinates. The kinematical system is 



r 






<t> 




v<f,/r 
dr + vl/r 


v r 




Vtf, 




a</> — VrV^/r 



(26) 



This new system (|26p makes use of the new acceleration formula in equation Q22p • Such a system can be solved using forward 
integration. Initial conditions for the photon have initial position at perihelion of 696000 kilometres about a mass with 
Schwarzschild radius of 2.95 kilometres, zero radial velocity, and tangential velocity v$ — c/n = Cy/Jl to the right. The speed 
of light is taken to be c = 300000 kilometres/second. Using Matlab's ODE45 routine (an explicit Runge-Kutta 4-5 method) 
produces the path shown in Fig. ((3} for the section shown; beyond this, the path is almost a straight line and so is not shown. 
The slope of the line between the last integration point before the photon reaches earth orbit, and the first point after is 
1.74851634 ± 10 -8 arcseconds. The time delay is calculated in the same integration, and is found to be 129.089609 ± 10 -6 
^seconds. The uncertainty is due to limitations on the precision of Matlab's ODE45 routine. Both of these values correspond 
well with the predictions from general relativity (as calculated above), more closely than the first order approximations, and 
use of a higher precision computation will allow a more accurate result, should such be required. Accuracy beyond first order is 
not commonly required, but these results give confidence that the acceleration vector presented here does accurately embody 
the effect on the photon due to a single Schwarzschild- type mass. 

Using the values found for position and radial and tangential velocities when the photon reaches earth orbit, we can send 
the photon along the path from earth orbit past the sun and back out to earth orbit. The central section is shown in Fig. ((4]). 
Note that the scales differ by a factor of 10 6 . 

The values for the deflection and the Shapiro delay calculated by this forward integration ray shooting are compared 
with the predictions from general relativity and the usual first order approximations(as calculated earlier) in Table [T] and 
demonstrate that the kinematic ray shooting method presented in this paper is an accurate representation of the effect of the 
gravitational field of a single Schwarzschild body on the motion of a photon, giving us some confidence in using this method 
in approximating more complicated systems. 
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Figure 3. Path of photon from perihelion, plotted in Matlab using 2d kinematic ray shooting 




0.5 1 1.5 2 2.5 3 

x (kilometres from perihelion) 6 




4.2 Magnification Maps 

When considering a multi-body system, such as a planetary system, it must be stressed that there is no known metric. That is, 
there is no known exact solution to Einstein's equations for such a system. Some sort of approximation is therefore required. 
Use of the 'weak-field metric' is one such approximation, as is the addition of the deflection angles due to each body (the 
method generally used in microlensing models). Here, we choose to approximate by adding the acceleration components due 
to each body in the system. 

Having tested the radial and tangential acceleration components described above, it is a simple matter to set up a three 
dimensional ray tracing system for a planetary system. At each integration point along each ray, the acceleration components 
due to each massive body are calculated. This is done by a translation to put the massive body at the origin, followed by three 
rotations to place the photon's position vector and its velocity vector in the same plane as the massive body, with 6 — it/2. 
The radial and tangential velocities are then used to calculate the radial and tangential acceleration components. The three 
rotations are then reversed, and the resultant cartesian acceleration components are added to the acceleration components 
due to any other masses in the system to determine the overall acceleration of the photon. 

As a very simple (and artificial) example of this process, we first consider a two dimensional system, consisting of two very 
massive bodies in close proximity, and plot the path of the photon through this binary system. The sma ller and larger stars 



repres ent bodies of 20 million and 50 million solar masses respectively (similar to the system described in lBoroson and Lauer 



(2009), although in the present example we imagine that the system has decayed to the point where the black holes are only 
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Figure 5. Path of a photon through a supermassive binary black hole system. The black holes arc indicated on the diagram with 
asterisks, and are located at (-946.08,-315.36) and (0,23.652). 
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a billion kilometres apart). This example is purely to demonstrate the versatility of the present approach. We are making the 
gross simplification that the black holes are stationary throughout the period when the ray is passing through the system, 
and so the acceleration of the two masses towards each other is therefore being ignored. For the purpose of demonstrating the 
procedure used here, we are ignoring such limitations. The present model is clearly a coarse approximation in this extreme 
case. The light ray comes in from the far left, and is deflected by the summed acceleration components due to each mass. The 
path is shown in Fig. J5). 

We now consider a more commonly modeled lensing system, the planetary lens. In describing a lensing system, it is 
common to use a parameter called the Einstein radius. This is the angular radius at which observers perfectly aligned with 
the point source and point lens would see a ring of light about the lens. Specifically, the system is designed as follows: a point 
source is at (—8000, 0, 0) and the observer is at (+8000, 0, 0). The lens star is placed at the origin, having Schwarzschild radius: 
r s = 99 * 10 -8 . An orbiting planet is placed at (0,0.1208, 0) (in microlensing terms, it is at 1.35 times the Einstein radius), 
with r s = 1 * 10~ 8 . For simplicity in this model, we ignore the motion of the planet. Rays are sent through the system, near 
the Einstein radius, and in the vicinity of the planet. Due to symmetry in the cases here, it is only necessary to calculate 
half the rays and plot the result both above and below the axis of symmetry. During the numerical integration, each ray is 
broken into several small sections, with the size of each section becoming smaller as the photon nears the lens star. This is 
important to ensure that the integration routine does not take too large a step and miss the strong deflection altogether. 
For the simulations presented in this paper, each ray is broken into 36 segments. The result of this process is a light density 
map, or magnification map (that is points in (y, z) where the rays cross the plane x = 8000). These simulations were run on 
MatLab version 6.1, under Windows XP, on an Athlon 64 3500+ processor with 1GB of ram. Running times are given as an 
indication, but no measures have been made to optimise the code for efficiency. 

Fig. (j6]) shows the magnification map produced when a rectangular array of rays (222 by 205) is sent through the lens 
system. The bending of the light towards the planet is clearly visible. This, combined with the bending cause d by the lens 
star, p roduces the characteristic diamond shape for a system with a single planet outside the Einstein radius jWambsgansi ' 



w) 



(1997)). The running time was 26 hours. In Fig. (:7;). many more rays have been used (approximately 106.000 rays). In order 
to view the resulting density, it is necessary to col our or shade r egions according to how many rays pass through each small 
area. The code used to do this was 

l smoothhist2d',[Perkin] diooi)- The running time was approximately 50 hours. This result 
clearly shows the caustic diamond structure. 

This diamond structure is the expected shape for the magnification map, and suggests that the method used here can be 
considered as an alternative method for modeling a thin lens. 
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Figure 6. Caustic structure due to planet with mass 1% of star's mass, located at 1.35 Einstein radii, using over 40,000 rays 
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Figure 7. magnification density plot; same parameters as in Fig. JB), using 106000 rays 
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5 CONCLUSION AND DISCUSSION 

We have considered the path of a photon near a Schwarzschild-type body. Using the Schwarzschild metric, a new refractive 
index has been derived. Integrating the angle along the path gives the total deflection angle, and integrating the new index 
along the path gives the travel time delay. These values for the delay time and deflection can be calculated, using the formulae 
here, to arbitrary precision. This is because these formulae are derived directly from the Schwarzschild metric. As a check, 
it was shown that the standard first order approximations used for the deflection and delay agree with these results to first 
order. 

A new formula for the acceleration of a photon was derived by combining the Schwarzschild metric with the path equation, 



A kinematical approach to gravitational lensing 11 



and differentiating. This new acceleration formula was tested with a ray shooting approach, using the new refractive index 
to provide the initial conditions for the velocity components. The deflection and delay values were found to be in excellent 
agreement with the precise values calculated earlier. 

By making the approximation that the acceleration on the photon is the sum of the individual acceleration components 
due to each massive body, a simple microlensing model was developed to demonstrate a use of the new acceleration formula 
for a binary system. Sample light fields on an observer plane have been computed using this new approach, and reproduce 
the expected figures. No attempt has been made here to speed computations, since that was not the purpose of the present 
work, but future developments may address such issues of computational efficiency. 

In summary, this work provides a 'classical' way of accurately describing the gravitational effect on a photon due to 
a single mass, and provides an alternative method for approximating the course of a photon through a complicated mass 
distribution. The authors believe that the approach presented here provides an insight into the effect of gravitating bodies 
on light rays that can be grasped without requiring a deep understanding of general relativity, and yet is still quantitatively 
accurate for a single mass, and can be used for approximating more complicated systems. As for the standard thin lens 
'deflection angle' method, this approach to gravitational lensing may be used by applied mathematicians, computer modellers 
and others without requiring specialist knowledge of general relativity. Because this approach retains the 'delay' information 
as well as the deflection, it might conceivably be useful in analysing systems where the time delay plays a role, such as a pulsar 
source being lensed, should we observe such an event. It is also hoped that the formulae presented here will prove useful in 
producing models of more complicated mass distributions, such as galaxy clusters. Such models could be produced using the 
same method used here, simply by adding more bodies to the model. 

ACKNOWLEDGEMENT 

The authors are indebted to O. Wucknitz for constructive comments on an earlier draft of this paper. 



REFERENCES 

Boroson T.A., Lauer T.R., 2009, Nature Lett. Vol 458, 53-55. 

Capozziello S., de Ritis R., Man'ko V.I., Marino A. A., Marmo G., 1997, Physica Scripta Vol 56, 212-220. 
Carroll S.M., 2004, Spacetime and Geometry. Addison Wesley, San Francisco 
Ellis R., 2010, Phil. Trans. R. Soc. A, 368, 967-987 

Kenyon I.R., 1990, General Relativity. Oxford University Press, Oxford 

Mao S., Kerins E., Rattenbury N.J., 2007, Proceedings IAU Symposium, No. 249 

Misner C. W., Thorne K. S., Wheeler J. A., 1973, Gravitation. Freeman, San Francisco, Chapter 23 

Newton I., 1704, Opticks, Fourth Ed., Reprint 1931, G. Bell and Sons, London 

Paczynski B., 1986, ApJ, 304, 1-5 

Perkins P., 2006, MatLab Central File Exchange, 

http://www.mathworks.com.au/matlabcentral/fileexchange/13352-smoothhist2d 
Schneider P., Ehlers J., Falco E., 1992, Gravitational Lenses, Springer- Verlag, New York 
Shapiro I.I., 1964,Physical Review Letters 13 (26), 789-791 
Sumi T., Bennett D.P., Bond LA. et al, 2010, ApJ, 710:16411653 
Wambsganss J., 1997, Mon. Not. R. Astron. Soc, 284, 172-188 
Weinberg, S., 1972, John Wiley & Sons Inc. 
Zabel S. A., Peterson J. B., 2003, ApJ, 594, 456-463 



